(Vector and) Raster Data in R

Why care about data types and formats?

There are differences in how spatial information is stored, processed, and visually represented.

• Different commands for data import and manipulation
• Spatial linking techniques and analyses partly determined by data format
• Visualization of data can vary

So, always know what kind of data you are dealing with!

Vector and raster data

Sources: OpenStreetMap / GEOFABRIK (2018), City of Cologne (2014), and the Statistical Offices of the Federation and the Länder (2016) / Jünger, 2019

Representing the world in vectors

Vector data: geometries

Vector data: attribute tables

Only geometries means that we do not have any other information.

We must assign attributes to each geometry to hold additional information \(\rightarrow\) data tables called attribute tables.

• Each row represents a geometric object, which we can also call observation, feature, or case
• Each column holds an attribute or, in ‘our’ language, a variable

Vector data: attribute tables

File formats/extensions

• GeoPackage .gpkg
• Shapefile .shp
• GeoJSON .geojson
• …
• Sometimes, vector data come even in a text format, such as CSV

Welcome to simple features

Load a vector data file

The first step is, of course, loading the data. We want to import the .geojson file for the administrative borders of the whole world called World_Countries.geojson.

# load library
library(sf)

# load data
# source: https://hub.arcgis.com/datasets/esri::world-countries-generalized/
hello_world <- sf::read_sf("./data/World_Countries.geojson")

What is this thing?

hello_world

Simple feature collection with 251 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89 xmax: 180 ymax: 83.6236
Geodetic CRS:  WGS 84
# A tibble: 251 × 6
     FID COUNTRY             ISO   COUNTRYAFF  AFF_ISO                  geometry
   <int> <chr>               <chr> <chr>       <chr>          <MULTIPOLYGON [°]>
 1     1 Afghanistan         AF    "Afghanist… "AF"    (((61.27655 35.60725, 61…
 2     2 Albania             AL    "Albania"   "AL"    (((19.57083 41.68527, 19…
 3     3 Algeria             DZ    "Algeria"   "DZ"    (((4.603354 36.88791, 4.…
 4     4 American Samoa      AS    "United St… "US"    (((-170.7439 -14.37555, …
 5     5 Andorra             AD    "Andorra"   "AD"    (((1.445836 42.60194, 1.…
 6     6 Angola              AO    "Angola"    "AO"    (((23.47611 -17.62584, 2…
 7     7 Anguilla            AI    "United Ki… "GB"    (((-63.16778 18.16445, -…
 8     8 Antarctica          AQ    " "         " "     (((-180 -84.30535, -179.…
 9     9 Antigua and Barbuda AG    "Antigua a… "AG"    (((-61.73806 16.98972, -…
10    10 Argentina           AR    "Argentina" "AR"    (((-71.85916 -41.01128, …
# ℹ 241 more rows

We can already plot it

plot(hello_world["FID"])

This is the bounding box

Inspect your data: classics

There are no huge differences between the file we just imported and a regular data table.

# head of data table
head(hello_world, 2)

Simple feature collection with 2 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 19.28854 ymin: 29.40611 xmax: 74.91574 ymax: 42.66035
Geodetic CRS:  WGS 84
# A tibble: 2 × 6
    FID COUNTRY     ISO   COUNTRYAFF  AFF_ISO                           geometry
  <int> <chr>       <chr> <chr>       <chr>                   <MULTIPOLYGON [°]>
1     1 Afghanistan AF    Afghanistan AF      (((61.27655 35.60725, 61.29638 35…
2     2 Albania     AL    Albania     AL      (((19.57083 41.68527, 19.58195 41…

Inspect your data: spatial features

Besides our general data inspection, we may also want to check the spatial features of our import. This check includes the geometric type (points? lines? polygons?) and the coordinate reference system.

# type of geometry
sf::st_geometry(hello_world) 

Geometry set for 251 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89 xmax: 180 ymax: 83.6236
Geodetic CRS:  WGS 84
First 5 geometries:

Inspect your data: spatial features

Each polygon is defined by several connected points to build an enclosed area. Several polygons in one data frame have the sf type multipolygons. Just as the world consists of several states, the polygon of the world consists of several smaller polygons.

# extract the simple features column and further inspecting 
attr(hello_world, "sf_column") |> 
  dplyr::pull(hello_world, var = _) |> 
  dplyr::glimpse()

sfc_MULTIPOLYGON of length 251; first list element: List of 1
 $ :List of 1
  ..$ : num [1:708, 1:2] 61.3 61.3 61.4 61.4 61.4 ...
 - attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"

Inspect your data: spatial features

Remember: The Coordinate Reference System is critical. A crucial step is to check the CRS of your geospatial data.

# coordinate reference system
sf::st_crs(hello_world) 

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Spatial joins

Often we want to combine geospatial datasets from different sources. The choice of an appropriate method depends, once again, on the data format. In the world of vector data, we usually speak of “spatial joins” when linking two vector datasets. Let’s pull in more spatially narrowed data.

german_districts <- 
  sf::st_read(
    "./data/VG250_KRS.shp", 
    as_tibble = TRUE
  )

plot(german_districts["AGS"])

Spatial joins

We can craft a toy example to demonstrate spatial joins using the synthetic survey geocoordinates. This is how they look.

points <- 
  readRDS(
    "./data/synthetic_survey_coordinates.rds"
  ) |> 
  sf::st_transform(sf::st_crs(german_districts))

plot(points)

Spatial joins

It would be nice, if we had regional identifiers in our points data but there aren’t any. This issue shows where spatial joins come in handy. We can use sf::st_join() to extract this information.

points_ags <-  sf::st_join(points, german_districts["AGS"])

points_ags

Simple feature collection with 2000 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 294988.7 ymin: 5261146 xmax: 904079.4 ymax: 6080921
Projected CRS: ETRS89 / UTM zone 32N
# A tibble: 2,000 × 3
      id           geometry AGS  
 * <int>        <POINT [m]> <chr>
 1     1 (853435.9 5892079) 12073
 2     2 (520533.4 5879316) 03357
 3     3 (430178.6 5689132) 05958
 4     4 (521463.3 5884327) 03357
 5     5 (525072.6 5397485) 08116
 6     6   (765168 5904791) 13071
 7     7 (475594.9 5803719) 05770
 8     8 (670010.2 5473386) 09574
 9     9 (669450.7 5362344) 09174
10    10 (827147.8 5388243) 09262
# ℹ 1,990 more rows

There’s more

In this course, we only use a subset of the options for spatial joins in sf. For more please refer to the official sf cheat sheet.

Finally: Transforming the CRS

Depending on our data sources, we have to change our CRS fairly often. That’s easily doable in sf using the sf::st_transform(). Here’s a closer look.

sf::st_crs(german_districts)$epsg

[1] 25832

Let’s pretend we want to use the EPSG code 4326–no problem!

german_districts_4326 <-  sf::st_transform(german_districts, 4326)

sf::st_crs(german_districts_4326)$epsg

[1] 4326

Difference to vector data

Data Structure:

• Other data format(s), different file extensions
• Geometries do not differ within one dataset

Implications:

• Other geospatial operations possible

Benefits:

• Can be way more efficient and straightforward to process
• It’s like working with simple tabular data

Visual difference between vector and raster data

What exactly are raster data?

• Hold information on (most of the time) evenly shaped grid cells
• Basically, a simple data table
• Each cell represents one observation

Metadata

• Information about geometries is globally stored
• They are the same for all observations
• Their location in space is defined by their cell location in the data table
• Without this information, raster data were simple image files

Important metadata

Raster Dimensions: number of columns, rows, and cells

Extent: similar to bounding box in vector data

Resolution: the size of each raster cell

Coordinate reference system: defines where on the earth’s surface the raster layer lies

Setting up a raster dataset is easy

input_data <- matrix(sample(1:100, 16), nrow = 4)

raster_layer <- terra::rast(input_data)

raster_layer

class       : SpatRaster 
dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 4, 0, 4  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
name        : lyr.1 
min value   :     1 
max value   :    95 

We can already plot it

terra::plot(raster_layer)

File formats/extensions

• GeoTIFF tif
• Gridded data .grd
• Network common data format .nc
• Esri grid .asc
• …
• Sometimes, raster data come even in a text format, such as CSV

Implementations in R

terra is the most commonly used package for raster data in R.

Some other developments, e.g., in the stars package, also implement an interface to simple features in sf. We will work with stars later in this course.

The terra package helps to employ zonal statistics. But the spatstat package is even more elaborated.

Loading raster tiffs (raw WorldPop data)

Loading raster data in terra is straightforward. For this purpose, we use the function terra::rast().

pop_grid_ger_2020 <- terra::rast("./data/deu_ppp_2020_1km_Aggregated.tif")

pop_grid_ger_2020

class       : SpatRaster 
dimensions  : 934, 1100, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.87375, 15.04042, 47.27458, 55.05792  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : deu_ppp_2020_1km_Aggregated.tif 
name        : deu_ppp_2020_1km_Aggregated 
min value   :                       0.000 
max value   :                    4513.521 

Loading raster tiffs (raw WorldPop data)

Let’s load another dataset for comparing structures, because that’s really easy in terra.

pop_grid_ger_2000 <- terra::rast("./data/deu_ppp_2000_1km_Aggregated.tif")

pop_grid_ger_2000

class       : SpatRaster 
dimensions  : 934, 1100, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.87375, 15.04042, 47.27458, 55.05792  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : deu_ppp_2000_1km_Aggregated.tif 
name        : deu_ppp_2000_1km_Aggregated 
min value   :                       0.000 
max value   :                    3367.343 

Compare layers by plotting

terra::plot(
  pop_grid_ger_2000, 
  range = c(0, 4000)
)

terra::plot(
  pop_grid_ger_2020, 
  range = c(0, 4000)
)

Simple statistics

Working with raster data is straightforward

• quite speedy
• yet not as comfortable as working with sf objects

For example, to calculate the mean, we can use the following:

terra::global(pop_grid_ger_2000, fun = "mean", na.rm = TRUE)

                                mean
deu_ppp_2000_1km_Aggregated 125.3752

Get all values as a vector

We can also extract the values of a raster directly as a vector:

all_raster_values <- terra::values(pop_grid_ger_2000)

mean(all_raster_values, na.rm = TRUE)

[1] 125.3752

Nevertheless, although raster data are simple data tables, working with them is a bit different compared to, e.g., simple features.

Combining raster layers to calculate new values

pop_grid_ger_2000_2020 <-  
  pop_grid_ger_2020 - pop_grid_ger_2000

pop_grid_ger_2000_2020

class       : SpatRaster 
dimensions  : 934, 1100, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.87375, 15.04042, 47.27458, 55.05792  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : deu_ppp_2020_1km_Aggregated 
name        : deu_ppp_2020_1km_Aggregated 
min value   :                   -901.1108 
max value   :                   1609.1487 

Loading raster data with stars

Later, we will work with multidimensional raster data cubes. While it is possible to work them a bit also in terra, stars is way more elaborated. Just that you already can see, how to import raster data with stars, here’s an example.

pop_grid_ger_2020_stars <- 
  stars::read_stars("./data/deu_ppp_2020_1km_Aggregated.tif")

pop_grid_ger_2020_stars

stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                                      Min.  1st Qu.   Median    Mean  3rd Qu.
deu_ppp_2020_1km_Aggregated...  0.01720133 7.404476 14.82611 60.1534 32.34949
                                    Max.  NA's
deu_ppp_2020_1km_Aggregated...  2277.353 83687
dimension(s):
  from   to offset     delta refsys point x/y
x    1 1100  5.874  0.008333 WGS 84 FALSE [x]
y    1  934  55.06 -0.008333 WGS 84 FALSE [y]

stars objects can also be plotted

plot(pop_grid_ger_2020_stars)

Finally: Transforming the CRS

As mentioned before, we may have to change the CRS also of our raster data. Going back to terra, this is not called ‘transforming’ but ‘projecting’ (which makes sense, right?).

terra::crs(pop_grid_ger_2020, describe = TRUE)$code

[1] "4326"

Let’s pretend we want to use the EPSG code 25832–no problem!

pop_grid_ger_2020_25832 <- terra::project(pop_grid_ger_2020, "EPSG:25832")

terra::crs(pop_grid_ger_2020_25832, describe = TRUE)$code

[1] "25832"